library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix

rm(list = ls())

# load data
load(file = 'chDataWomenDistrict.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]

# "Gonzalo Fuenzalida" and 'Javier Macaya' changed district after the reform
ch_reform <- subset(ch_reform, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTerm <- subset(ch_refTerm, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_pre_post <- subset(ch_pre_post, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChanges <- subset(ch_refTermMChanges, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
table(subset(ch_reform, bothMag == 1)$bill_author, subset(ch_reform, bothMag == 1)$reform)


################################# 
####### Models using M ##########
#################################

# M1: G*M C1,C2,C3 (all data)
data_m1 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman + n_of_women
	~ session, data = ch_refTerm)
data_m1 <- cbind(data_m1, bill_author = ch_refTerm$bill_author)
m1 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman  + n_of_women - 1, data = data_m1)
vcov1 <- cluster.vcov(m1, cluster = data_m1$bill_author)
s1  <- sqrt(diag(vcov1))

# M2: G*M C1,C2,C3 (M changes at the middle of C2)
data_m2 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman + n_of_women
	~ session + reform, data = ch_refTermMChanges)
data_m2 <- cbind(data_m2, bill_author = ch_refTermMChanges$bill_author)
m2 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman + n_of_women - 1, data = data_m2)
vcov2 <- cluster.vcov(m2, cluster = data_m2$bill_author)
s2  <- sqrt(diag(vcov2))


# M3: G*M C1,C3
data_m3 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman + n_of_women
	~ session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022')))
data_m3 <- cbind(data_m3, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))$bill_author)
m3 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman + n_of_women  - 1, data = data_m3)
vcov3 <- cluster.vcov(m3, cluster = data_m3$bill_author)
s3  <- sqrt(diag(vcov3))

# M5: G*M C1,C2,C3 (only repeated legislators)
data_m5 <- fixest::demean(women_pct + magAnalysis + magAnalysisWoman + n_of_women
	~ bill_author + session, data = subset(ch_refTerm, allSessions == 1))
data_m5 <- cbind(data_m5, bill_author = subset(ch_refTerm, allSessions == 1)$bill_author)
m5 <- lm(women_pct ~ magAnalysis + magAnalysisWoman + n_of_women  - 1, data = data_m5)
vcov5 <- cluster.vcov(m5, cluster = data_m5$bill_author)
s5  <- sqrt(diag(vcov5))

# M6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
data_m6 <- fixest::demean(women_pct + magAnalysis + magAnalysisWoman + n_of_women
	~ session + reform + bill_author, data = subset(ch_refTermMChanges, allSessions == 1))
data_m6 <- cbind(data_m6, bill_author = subset(ch_refTermMChanges, allSessions == 1)$bill_author)
m6 <- lm(women_pct ~ magAnalysis + magAnalysisWoman + n_of_women  - 1, data = data_m6)
vcov6 <- cluster.vcov(m6, cluster = data_m6$bill_author)
s6  <- sqrt(diag(vcov6))

# M7: G*M C1,C3 (only repeated legislators)
data_m7 <- fixest::demean(women_pct + magAnalysis + magAnalysisWoman + n_of_women
	~ bill_author + session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1))
data_m7 <- cbind(data_m7, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1)$bill_author)
m7 <- lm(women_pct ~ magAnalysis + magAnalysisWoman + n_of_women  - 1, data = data_m7)
vcov7 <- cluster.vcov(m7, cluster = data_m7$bill_author)
s7  <- sqrt(diag(vcov7))

# M8: C3
data_m8 <- subset(ch_refTerm, session == '2018-2022')
m8 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman + n_of_women, data = data_m8)
vcov8 <- vcovHC(m8, 'HC2')
s8  <- sqrt(diag(vcov8))

# Generate Table I.33
stargazer(m1, m2, m3, m5, m6, m7, m8, 
	se = list(s1, s2, s3, s5, s6, s7, s8), 
	covariate.labels = c('M', 'Woman', 'M x Woman', '# of Women in the District'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'Yes', 'No', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--2014-2019  Cámara de Diputadas y Diputados de Chile (Women in the District)",
	label = 't:womenDistrict'
)

### Substantive Effect
## Scenarios: Women, M = 2
W2 <- data.frame(mag = 2, female = 1, magAnalysisWoman = 2 * 1, n_women = 0)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = 3, female = 1, magAnalysisWoman = 3 * 1, n_women = 0)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = 2, female = 0, magAnalysisWoman = 2 * 0, n_women = 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = 3, female = 0, magAnalysisWoman = 3 * 0, n_women = 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)


## Sims
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims5 <- rmvnorm(n = 1000, mean = coef(m5), sigma = vcov5)
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)
sims7 <- rmvnorm(n = 1000, mean = coef(m7), sigma = vcov7)
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)

## Calculate Predicted Values
# Women
pred1_W <- apply((sims1 %*% t(Wdata[2,])) - (sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_W <- apply((sims2 %*% t(Wdata[2,])) - (sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_W <- apply((sims3 %*% t(Wdata[2,])) - (sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred5_W <- apply((sims5 %*% t(Wdata[2, -2])) - (sims5 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred6_W <- apply((sims6 %*% t(Wdata[2, -2])) - (sims6 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred7_W <- apply((sims7 %*% t(Wdata[2, -2])) - (sims7 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred8_W <- apply((sims8 %*% t(cbind(1, Wdata)[2,])) - (sims8 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))

# Men
pred1_M <- apply((sims1 %*% t(Mdata[2,])) - (sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_M <- apply((sims2 %*% t(Mdata[2,])) - (sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_M <- apply((sims3 %*% t(Mdata[2,])) - (sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred5_M <- apply((sims5 %*% t(Mdata[2, -2])) - (sims5 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred6_M <- apply((sims6 %*% t(Mdata[2, -2])) - (sims6 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred7_M <- apply((sims7 %*% t(Mdata[2, -2])) - (sims7 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred8_M <- apply((sims8 %*% t(cbind(1, Mdata)[2,])) - (sims8 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))


# Building blocks of table 2
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred5_W[2,], 2), round(pred6_W[2,], 2),
	round(pred7_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred5_W[3,], 2), ', ', round(pred5_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred7_W[3,], 2), ', ', round(pred7_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred5_M[2,], 2), round(pred6_M[2,], 2),
	round(pred7_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred5_M[3,], 2), ', ', round(pred5_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred7_M[3,], 2), ', ', round(pred7_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableI34 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(tableI34) <- paste0('Model ', c(1:3, 5:8))
stargazer(tableI34,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")

